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Abstract: This research presents a distributed and formula-based bilateration algorithm 
that can be used to provide initial set of locations. In this scheme each node uses 
distance estimates to anchors to solve a set of circle-circle intersection (CCI) problems, 
solved through a purely geometric formulation. The resulting CCIs are processed to pick 
those that cluster together and then take the average to produce an initial node location. 
The algorithm is compared in terms of accuracy and computational complexity with a 
Least-Squares localization algorithm, based on the Levenberg-Marquardt methodology. 
Results in accuracy vs. computational performance show that the bilateration algorithm is 
competitive compared with well known optimized localization algorithms. 

Keywords: distributed-localization; wireless sensor networks; Least Squares (LS); 
optimization; bilateration 
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1. Introduction 

Recent advances in microelectronics have led to the development of autonomous tiny devices called 
sensor nodes. Such devices, in spite of their physical limitations, contain the essential components 
of a computer, such as memory, I/O ports, sensors, and wireless transceivers which are typically 
battery-powered. Once deployed (randomly or not) over a certain area, sensor nodes have the ability 
to be in touch via wireless communications with neighboring nodes forming a wireless sensor network 
(WSN). The great advantage of using WSNs is that they can be applied in important areas such as disaster 
and relief, military affairs, medical care, environmental monitoring, target tracking, and so on [1-3]. 
However, most of WSN applications are based on local events. This means that each sensor node needs 
to detect and share local phenomenons with neighboring nodes, implying that the location of such events 
{i.e. , sensor locations) are crucial for the WSN application. In this way, sensor self-positioning represents 
the first startup process in most WSN projects. It is well known that using a GPS in each sensor node 
represents the primary solution to infer position estimates. However, this option is not suitable to be 
considered in all nodes if parameters like size, price, and energy-consumption in a sensor node are of 
concern [4]. In order to optimize such parameters, a good option consists of reducing to a small fraction 
of sensors with GPS, and the remainder sensors (i.e., unknown sensors), commonly above 90% of total 
deployed sensors, should use alternatives to estimate its own positions like radio-frequency transmissions 
or connectivity with neighboring sensors [5-7]. 

In order to provide position estimates many localization algorithms have been proposed, coming 
from different perspectives as described in [8,9]. Basically, localization algorithms can be categorized 
according to range -based vs. range-free methods, anchor-based vs. anchor-free models, and distributed 
vs. centralized processing [10,11]. Range-based methods consist of estimating node locations (using a 
localization algorithm) based on ranging information among sensor nodes. Range estimation between 
pairs of nodes is achieved using techniques of Time-of-Arrival (ToA), Receive- Signal- Strength (RSS), 
or Angle-of- Arrival (AoA) [12]. This approach has the disadvantage of requiring extra-hardware in 
each sensor board, increasing the cost per sensor. However, as far as is known, this approach provides 
the best cost-accuracy performance in localization algorithms. A less expensive but more inaccurate 
alternative consists of using just connectivity among sensor nodes to estimate node locations, called 
range-free [13]. On the other hand, if position estimates are obtained by considering absolute references 
(e.g., sensors with GPS or Anchors), the resulted position estimates (also with absolute positions) 
will be closely related to such reference positions, called an anchor-based model. By the contrary, if 
no reference positions are used to estimate locations, relative coordinates will be obtained, called an 
anchor-free model. 

One of the most interesting and relevant aspects in WSN localization is associated with the way 
to compute the location of sensor nodes. For example, if all pairwise distance measurements among 
sensor nodes are sent to a central node to compute position estimates, the localization algorithm becomes 
centralized. This kind of central processing has the advantage of global mapping, but it has basically 
two important disadvantages which demerit its use in many cases when robustness and saving-energy 
have high priority in a WSN [14]. Some important centralized schemes are the next. In [15] an iterative 
descent procedure (i.e., Gauss-Newton method) is used in a centralized way to solve the Non-Linear 
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Least-Square (NLLS) problem. Another interesting centralized scheme was proposed in [16] where 
the WSN localization problem is modeled as linear or semidefinite program (SDP), and a convex 
optimization is used to solve problem. 

In contrast, when each sensor node estimates its own location using available information of 
neighboring nodes (e.g., range, connectivity, location, etc.), the localization process becomes distributed. 
Distributed processing is much less energy consuming in WSNs than centralized processing because 
centralized schemes need to collect relevant information from all nodes in the network which implies 
re-transmissions in multi-hop environments. Also, distributed algorithms are tolerant to node failures 
due to node redundancy. Thus, basically a distributed algorithm allows robustness, saving-energy, and 
scalability [14,17,18], which overcomes the limitations imposed by the centralized approach. In [19], 
a robust least squares scheme (RLS) for multi-hop node localization is proposed. This approach 
reduces the effects of error propagation by introducing a regularization parameter in the covariance 
matrix. However, the computational cost to mitigate the adverse effects of error propagation is too 
high at energy-constrained nodes. Similarly, [20] proposes two weighted least squares techniques 
to gain robustness when a non-optimal propagation model is considered however they failed to 
introduce a covariance matrix in the localization process that can effectively decrease the computational 
complexity. On the other hand, the authors of [21] propose a Quality of Trilateration method (QoT) 
for node localization. This approach provides a quantitative evaluation of different geometric forms of 
trilaterations. However, it seems to be that the main idea of this methodology depends on the quality 
or resolution of geometric forms {i.e., like image processing) which is impractical to be implemented in 
resource-constrained devices with limited memory and processing capabilities {i.e., nodes). 

In this paper, we analyze a range-based bilateration algorithm (BL) that can be used in a distributed 
way to provide initial estimates for unknown sensors in a wireless sensor network (our analysis consider 
that each unknown sensor can determine its initial position communicating directly with several anchors). 
In this case, each node uses a set of two anchors and their respective ranges at a time to solve a set of 
circle intersection problems using a geometric formulation. The solutions from these geometric problems 
are processed to pick those that cluster around the location estimate and then take the average to produce 
an initial node location. Finally, we present a computational/accuracy comparison between the BL 
algorithm, based on closed-formulas, and a classical Least Squares (LS) approach for localization, based 
on the iterative Levenberg-Marquardt algorithm (LM). 

The outline of this paper is as follows. In Section 2 we examine a popular ranging technique for 
WSNs used in our simulations. In Section 3 we explore the localization problem from the Least Squares 
point of view. In Section 4 we analyze in detail the BL algorithm. In Sections 5 and 6 we evaluate the 
accuracy and computational-complexity performance respectively between the bilateration algorithm vs. 
LS schemes. Finally, we present our conclusions. 

2. Ranging Techniques 

This section presents a brief overview of an existing ranging technique used to estimate the true 
distance between two sensor nodes using power measurements, called Received Signal Strength (RSS). 
This technique is popular because sensor nodes do not require special hardware support to estimate 
distances. As a first approximation, considering the free space path loss model, the distance djj between 
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two sensors St and sj can be estimated by assuming that the power signal decreases in a way that is 
inversely proportional to the square of the distance [X/afj). However, in real environments the signal 
power is attenuated by a factor d~^p. The path-loss factor is closely related to geometrical and 
environmental factors, and it varies from 2 to 4 for practical situations [22]. In noiseless environments 
the power signal traveling from a sensor sj to a sensor Si can be measured according to the relation [23] 

where the path-loss factor (r\ p ) depends directly on the environmental conditions. Pq is the received 
power at the short reference distance of do = lm from the transmitter. Also, Pq can be computed by the 
Friis free space equation [24]. The log-distance path loss model 

P ( l J \dB)=P 0 -m p \og^- (2) 

do 

measures the average large-scale path loss between sensors Si and sj. The actual path-loss (in dB) is a 
normally distributed random variable: 

P { 1 ]) ~9t{pi ] \<5 2 SH ) (3) 

where Osh is given in dB and reflects the degradations on signal propagation due to reflection, refraction, 
diffraction, and shadowing. It can be seen that the linear measurements and distance estimates have a 
log-normal distribution with a multiplicative effect on the measurements. The noisy range measurement 
Rij can be obtained from Equations (2) and (3) as 

Rij = \0^hT (4) 

3. Least-Squares Multilateration Localization Algorithms 

In this section, we describe multilateration schemes that provide solutions to the Least-Square (LS) 
problem for location estimates using noisy ranging information derived from ToA or RSS ranging 
techniques. Consider a set of N wireless sensor nodes S = {s\,S2,--- ,sn}, randomly distributed over 
a 2-D region whose locations are unknown. We represent these unknown locations with vectors 
Z/ = [x;,yi] . Further, we assume the presence of a set A = {ai,a2, . . . of M reference or anchor 
nodes with known position q ; = [xy, yj\ T . Anchor nodes, a,-, are equipped with GPS or a similar scheme 
to self localize. Also, for practical situations M <t^N with M > 2. We develop our discussion assuming 
a 2-D scenario, but it can be easily generalized to the 3-D case. 

Moreover, we assume that any sensor can estimate pairwise ranges with its neighbors using 
time-of-arrival (ToA) or radio signal strength (RSS) techniques [24]. Denote the range estimate between 
the node Si and anchor aj as 

Rij = djj + e t j (5) 
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where djj is the true distance between aj and Sj, and e ;j represents the measurement error introduced 
by environmental noise, propagation distortion, and the ranging technique. Then the solution to the 
localization problem for a node Si consists of minimizing the sum of certain weighted error-distance 
function e w (-) as follows: 

M 

Pi = arg min £ e w ( 1 1 q 7 - - x 1 1 - Rij) = arg min^F (x) (6) 

x j=1 

where p ; = [xj,y/] represents the most likely position for the sensor Si that minimizes f , || • || represents 
the Euclidean norm, and e w (x) represents a function that provides a specific weight to the argument x 
{i.e., error distance). For example, the function e2(x) = (x) 2 , the LS formulation, is commonly used to 
solve Equation (6) due its tractability and efficiency in both mathematical and computational analysis. 
The LS problem can be solved either by closed-form solutions or by iterative methods. Next we describe 
both methodologies in detail. 

3.1. Closed-Form LS Multilateration 

Closed-Form methods have the advantage of fast time processing, which is useful for constrained 
devices {i.e., motes) where the energy conservation represents one of the major concerns. However, this 
approach is also subject to inaccurate estimates due to noisy ranging measurements, so in most cases this 
approach is not a suitable option in real WSN scenarios where current ranging techniques are not able to 
provide the required accuracy on the ranging measurements. For example, Spherical Intersection (SX), 
Spherical Interpolation (SI), and Global Spherical Least Squares (GSLS) [25] can solve a nonlinear set 
of equations using closed-formulas. These approaches provide good accuracy in the estimated positions 
under conditions like small biases and small standard deviations, but they also provide meaningless 
estimates under noisy environments [26]. A more robust closed-form scheme consists of using the 
classical LS multilateration discussed next [19,27,28]. 

Consider that a sensor Si with Cartesian position p ; = [jc ; ,y ; ] r has already estimated its range Rjj to M 
anchors. For each anchor aj with position q 7 = [xj,yj] T , an equation ||q 7 - — p ; -|| 2 = R 2 - is generated as 
shown the next formulas: 

||qi -p ; || 2 = 4 {xi -xi) 2 + (yi -yi) 2 =R 2 a 

||q 2 -p ( i| 2 =4 (x2-x l ) 2 + (y 2 -y l ) 2 =R 2 2 

llqM-p/ll 2 =R 2 iM (xM-Xi) 2 + (y M -yi) 2 =R 2 M 

The system of Equations (7) can be linearized by subtracting the first equation (j = 1) from the last 
M — 1 equations arriving to a linear system that can be represented in a matrix form as 



Ap, = b 



(8) 
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where 



X2-x\ yi-y\ 
x3-x\ y-i-y\ 

xm-x\ ynt-yi 



(9) 



(M-l)x2 



Rl-Rl+x\+y\-x\-yl 

RiM-Rl+xj+yj-xl-yi 



(10) 



(M-l)xl 



Now the least squares solution to Equation (8) is to determine an estimate for p ; that minimizes 



/i P( )=min<ji||A P; -b 112 
Pi I 2 



min<J ^(Ap^-b) 7 (Ap,--b) 
Pi I 2 



(ID 



After some manipulations we obtain the following 

1 



1. 



/( Pi ) = min -p A' Ap, - p A 1 b + -W b 
Pi I 2 2 



and the gradient of / at p, is 

V/(p ; )=A r A Pi -A r b = 0 
which provides the estimate {i.e., normal equations) to Equation (8): 

p, = (A T A) _1 A r b 



(12) 



(13) 



(14) 



Solving for Equation (14) may not work properly if A r A is close singular, so a recommended 
approach is to use a Tikhonov regularization as follows: 
For > 0 (e.g., close to zero) 



fit (Pi) = min | ^ 1 1 Ap,- - b 1 1 2 + 1 1 1 p, 1 

= mm J -pf A r A Pi - pf A r b + ^b r b + ^pf Pi 



(15) 



Then the gradient of at p, is 



V//i(Pi) = A r Ap ; - A r b +m = 0 



(16) 



Sensors 2012, 12 



845 



Factorizing we arrive to a robust estimate for the LS problem where the idea is to modify eigenvalues 
to avoid working with zero eigenvalues [19,29]. 

Pi=(A T A+fH)- 1 A T b (17) 

3.2. Iterative LS Algorithms 

Iterative methods are usually employed either when large-data set of information need to be processed 
or when an exact solution to a certain problem is not feasible (e.g., non-linear systems of equations) [30]. 
Optimization techniques represent a good alternative to solve such non-linear equations using an 
iterative procedure. Optimization algorithms that solve Non-Linear Least-Square (NLLS) problems 
{i.e. , the WSN localization problem) have been extensively proposed where the Newton or Quasi-Newton 
methods are iteratively used to minimizing some residuals [15,31,32]. The next paragraphs describe two 
well known iterative algorithms that are used to solve the NLLS problem: the Levenberg-Marquardt 
(LM) and the Trust-Region-Reflective (TRR). 

Assuming that a node denoted s,, with Cartesian position p ; = [jc,-,y,] r , estimates its distance Rij to 
M anchors denoted a.j, with positions q ; = [xj : yj] T , with j = 1, . . ., M. Consider the following residual 
error vector: 

/?;i-||p,--qi|| 



R(pO 

RiM — 1 1 P/ — Qm| 

Therefore, to find the more likely position of p ; , the program 

>r, 



(18) 



min/(p f )=min -R( Pi ) y R( Pi ) (19) 
P; Pi \2 J 

is solved, which is the least squares problem. 

To solve Equation (19) we employ the TRR algorithm and the LM algorithm. The TRR algorithm 
uses a sub-space trust-region method to minimize a function f(x). Here, approximations to / inside 
of a trust-region are iteratively required. The three main concerns in this algorithm are how to choose 
and compute the approximation to the function, how to choose and modify the trust region, and, finally, 
how to minimize over the sub-space trust-region. Even though the TRR algorithm provides an accurate 
solution for the WSN initial estimates, it is expensive (computationally speaking) for constrained sensor 
nodes [33]. 

On the other hand, the LM algorithm uses the search direction approach (a mix between the 
Gauss-Newton direction and the steepest descent direction) to find the solution to Equation (19). This 
algorithm outperforms the simple gradient descent methodology [34], and also it avoids dangerous 
operations with singular matrices as the pure Newton method does, so this methodology represents 
a good algorithm for comparison with the bilateration approach due to its robustness, speed, and 
accuracy [35]. Following the procedure presented in [29], Equation (19) can be solved by the Line 
Search Levenberg-Marquardt methodology shown in Algorithm 1, where ||-|| is the i-2 norm, I is the 
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identity matrix, is the estimated distance between the mote Sj and the anchor aj, J(p&) represents the 
Jacobian of R(p&) at the iteration k, and M/(pfc) is the merit function given by 

^R T (p,)R(p,) (20) 
The derivative of the merit function at the iteration k is 

M , / (p Jt )=J r (pjfc)R(pjfc) (21) 
Aim is the Levenberg-Marquardt direction, 

W = p||J r (p*)R(p*)|| (22) 



where p G (0, 1), and finally 



i M 

po = ^Lq. (23) 



provides the initial guess required for the TRR and LM iterative algorithms. 



Algorithm 1 Levenberg-Marquardt methodology. 

Require: an initial position po 
Ensure: a solution p^+i 
1: Initialize: k = 0,x = Threshold, p = 0.05 

2: do 

3: Solve: (J r (p fc )J(p^) +^I)A LM = -J r (Pfc)R(p^) 

4: Find the sufficient decrease (Armijo's condition): 

5: such that a k = {\)'for t = 0, 1,... 

6: satisfies M f (p k + a k A LM ) < M / (p yt ) + 1(T V^/P^Alm 

7: Update position: p^ + i = p + oc^Alm 

8: Update W : = P || J 2 " (P/t)R(P/t) || 
9: while(||p fe+ i -p^li < xor k< 100) 



4. A Bilateration Localization Method 

In this section we present the bilateration method for WSN localization which can be used as the 
initialization step for iterative localization schemes. This algorithm avoids iterative procedures, gradient 
calculations, and matrix operations that increase the internal processing in a constrained device. This 
research was done independently of the work presented in [36]. Even though both schemes share 
the same idea (i.e., bilateration), the procedure and the scope of both works are different as shown in 
Subsection 4.1. We show that it is possible to obtain a position estimate by solving a set of bilateration 
problems between a sensor node and its neighboring anchors, and then fusing the solutions according to 
the geometrical relationships among the nodes. Our aim is to find a scheme that can be deployed on a 
computationally constrained node. We argue that bilateration is an attractive option as the localization 
problem is divided on smaller sub-problems which can be efficiently solved on a mote. Next we 
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start our development by introducing the typical assumptions and definitions considered in a WSN 
localization problem. 

Let us define anchor subsets C A such that A^ = {a 7 -,a^} with j ^ k. Hence, there is a total of 
Q = (^f) anchor subsets. Without loss of generality, consider the case for one node Sj that receives from 
a subset Ayfc the anchor positions q 7 and q^, and computes the respective ranges Rtj and 7?,^ using RSS 
or ToA measurements. A possible geometrical scenario for this configuration is shown in Figure 1. We 
can appreciate from this example that the range estimates Rtj is larger than d\j and is shorter than d^. 
Now, consider the two range circles shown in the figure; one with its origin at qj and radius Ry, and the 
second with center in and radius R%, Next, define the two circle intersection points as gj and gj , 
where gj k is the reflection of gj k with respect to the (imaginary) line that connects q 7 and q^. In this case, 
the superscript jk represents the anchor subset A^. To simplify our discussion, we drop the superscripts, 
and only use them when more than one anchor subset is involved in our discussion. 

Figure 1. Sensor si finding its two feasible solutions (g ; ,g ; ) based on the anchors locations 
(Ok, q ; ) and their respective anchor range measurements fi?^,/?^) • 




...ir' 



In our approach, node Si determines the circle-circle intersections (CCI) g/ and g ; - by solving the 
closed-form expression reported in [37]. For instance, consider the two right triangles formed by the 
coordinates (q/,g ; ,f f ) and (qfe,g/,f?) in Figure 1, which satisfy the following relationships: 

dl + h 2 =Rjj (24a) 
d 2 kt + h 2 = Rl (24b) 

respectively. The distance dj t can be obtained by solving for h 2 in Equations (24a) and (24b): 

4-4=4-4 (25) 

and letting d = ||q ; - — q,t|| = dp +dk t resulting in 



R 2 j-d 2 =R 2 k -(d-d Jt ) 2 (26a) 
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Ri = R 



ij 



ik 



d 



d z + 2-d-d 

2 
ik 



Jt 



Jt 



Rfj-Rl + d 2 



2-d 



where the position f t = [x tl y t ] T is obtained as follows 



(26b) 
(26c) 



Finally, the circle intersection g,; = [x,-,y,] r is computed as 



(27) 



xi=x t ±2^-yj) 
yi=ytTj{xk-Xj) 



(28a) 
(28b) 



where q 7 = [xj, yj] T , (\k = [xfayk] > an d h i s easily obtained from Equation (24). The complementary 
signs of Equations (28a) and (28b) are used to obtain the solution for g, . 

Each node Si applies the CCI procedure using all Q subsets A^. For instance, gj k and gj k are obtained 
from the subset Aj^, gf and gj are obtained from the subset A«, and so on. Hence, a sensor node will 
have 2Q possible initial position estimates where half are considered mirror solutions which should be 
eliminated through the selection process described next. Geometrically, we expect that the true location 
will be located around the region where solutions form a cluster (i.e., half of the circle intersections 
should ideally intersect at the solution). Let us to consider the example shown in Figure 2. 

Figure 2. Sensor Si getting its initial estimation from three non-collinear anchors 

(aj,ak,a£). 




There are three anchors named aj, and ai and a node si that needs to be localized. The range 
estimate Rij is larger than du, the range estimate R^ is shorter than d^, and the range estimate Rj£ is 
shorter than du. Hence, Sj computes a set of of six location candidates given by {g J j k ,gj k , gj ,g{ , g^, gf^}- 
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As seen in the figure, all the mirror circle intersection estimates will tend to be isolated while the correct 
circle intersections will tend to cluster around the node location. For example, to decide between gj and 
gj candidate positions, generated using the anchors (a^a^), the sensor Sj obtains the minimum Square 
Euclidean sum from the location gj to each pair of candidate positions as follows: 



\\f = mm 



jk ji 

gi -sj 



jk -ji 

g; -g/ 



+ min 



gf-gf 



,j. k -o¥ 



(29) 



On the other hand, the sensor Si also obtains the minimum Square Euclidean sum from the location 
gj to each pair of candidate positions as follows: 



cp = mm 



-jk ji 

8/ Si 



-jk -ji 

g; -gj 



+ min 



if-gf 



5; fe! 



(30) 



jk — jk 

Finally, the lowest value of \|/ and cp helps to decide between choosing gj or gj . The process is 
repeated for all Q solution pairs to generate a set of disambiguated locations. 



Algorithm 2 General code used by every sensor Sj to get its initial position estimate p". 

Require: qk,Rik,with{k<r- l,...,M},and<2«- (^) 
Ensure: 
1 

2 
3 
4 

5 

6: 
7: 



0 1 



9: 
10: 

11: 

12: 

13 
14 
15 
16 
17 
18 
19 
20 
21 



Initialize: T <— [0,0] r 

for each subset A^ G ( 2 ) two-anchor subsets do 
y^0 
(p^0 

(gf^gf^ CC/ (qj,qk,Rij,Rik) {Return the two circle intersections} 
for each subset Ag m ^ Ajk G ( 2 ) two-anchor subsets do 

(gf",gf n ) «— CCI(q£,q m ,Ra,Ri m ) {Return the two circle intersections} 



jk 

Si - U 

jk 

sj -g 



Vl <- 

V2 

\|/ ^— \j/+min (vi , V2) {Return the minimum between vi and V2} 

jk „lm 



W2 



Si -Si 

-jk 

gj "g 



(111 



cp <— cp+min {w\,W2) {Return the minimum between w\ and W2} 
end for 

if (\|/ < cp) then 

T^T + gf 
else 

T^-T + gf 
end if 
end for 

p i ^ Q 

1 The algorithm omits the special cases where there are no circle intersections. The procedure for these instances 
is described further in the text. 
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Referring to our example, once node Sj removes the mirror locations, then an estimate of the node 
position can be formed by taking the average of the disambiguated set G = {g;f\g/^,gf }• 

The complete bilateration scheme is described in Algorithm 2. This is a distributed localization 
algorithm in the sense that each node can implement it and determine its position estimate, given the 
anchor positions and the range estimates Rij between each node and all the anchors. Since Algorithm 2 
uses only anchor measurement, it can be used as an initialization step to generate a set of position 
estimates that can be used with algorithms that integrate more information from anchor and non-anchor 
nodes (i.e., iterative distributed algorithms). 

Figure 3. Possible cases where the triangle inequality is not satisfied, (a) Case A; (b) Case B. 















i 








' or/ 








/ \ "A* 

















< 

q 


\ 



..••XX-. 



(a) 



(b) 



There are some anomalous cases which should be considered in the bilateration algorithm. In order 
to get its initial estimation p9, it is essential that every sensor Si gets the two location estimations from 
each one of the Q subsets even if the solutions are not feasible. For example, assume the two special 
cases shown in Figure 3. If we consider the left-side case on the figure, Rn is shorter than dik, and R^ is 
shorter than du, clearly the triangle inequalities are not satisfied since 



Rij +Rik 




> 


ll^-^ll 




Rij+H 


-q*|| 


> 


Rik 


(31) 




-q*|| 


> 


Ru 





As a consequence, the sensor Si will not be able to find any solution to Equation (28). In other words, if 
the two circles do not intersect with each other, it will not be feasible to find the circle intersections g, and 
g ; . Therefore, a relaxed estimation should be generated as described next. Considering that \\qj — q^|| is 
a constant distance between the anchors in set A^, the node Si takes two steps to estimate the locations 
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gj and gj . First, a location xi is obtained by fixing and making Rjj = | q ; - — — in order to 
satisfy the triangle inequality. Next, the sensor s; should use the CC7 procedure to solve for xi . Similarly, 
a second location estimate X2 is obtained by fixing Ry, choosing Rjk = \ \\c[j — qk\ \ —Rij\ to satisfy the 
triangle inequality and solving the problem through the CCI procedure. Finally, both gj k and gj are 
generated as the average X '^ X2 implying that when the triangle inequality is not satisfied, there will be 
a single solution that falls over the line y. A similar procedure can be derived for the second case as 
depicted in Figure 3. 

4.1. Comparison with the Previous Bilateration Scheme 

As described before, the research reported in [36] is focused on a distributed bilateration scheme that 
finds initial estimates. Using two anchors at a time, each sensor node s, finds two possible candidates 
{i.e., circle intersections). If sufficient anchors are available, the sensor node si averages the cloud 
of candidates which tend to be close to each other. The average of such candidates provides the 
initial estimate. 

As can be seen the general idea for this approach is quite similar to our bilateration approach. 
However, there are relevant differences between the two schemes that should be taken into account. 
These differences make that our bilateration approach be an alternative for the scheme proposed in [36]. 
For example, one of the differences is that [36] does not take into account special cases when a sensor 
Si is not able to compute circle intersections of two anchors {i.e., the circles are not in touch) as shown 
Figure 3. Therefore, under this perspective this scheme is limited to naive scenarios in which estimated 
distances between sensors and anchors should have good accuracy. Thus, noisy RSS measurements, 
commonly used in realistic scenarios, may not provide useful information for this scheme. Hence, if a 
sensor si is not able to find sufficient circle intersections from two neighboring pair of anchors at time, 
the localization process will fail. In our case, the bilateration scheme is able to obtain initial estimates 
under the most severe scenarios {i.e., not circle intersections). 

Another important aspect to consider in [36], is the use of a threshold, 8, which reduces the number 
of possible candidate positions, making this approach more selective. However, the value of 8 is hard to 
determine in practice and also it does not guarantee good results in noisy environments. Additionally, 
in [36] each sensor node Si should create a table of its neighboring anchors. All anchors have a specific 
position inside of the table, and they are weighted by the sensor st according to the candidate positions 
that they generate. The value of 8 is used to select a certain group of candidate positions. The anchors 
are weighted according to the candidates that they generated. Finally, all tables are broadcast by sensors. 
Once all sensors have received the anchor tables of their neighbors, they run a post-processing stage to 
determine which anchors are more reliable than others. These anchors are used to obtain initial estimates. 
As can be seen, the drawback of this approach are extra wireless transmissions required to share anchor 
tables among sensors. In our case we present an extension of the earlier BL algorithm which avoids any 
kind of wireless transmission with the goal to save energy. Finally, we should remark that we are using 
a sorting algorithm (lines 10,13, and 15-18 of the Algorithm 2) to determine initial positions. Analysis 
results shown in next section demonstrate that the alternative BL algorithm is competitive in comparison 
with well known accurate and efficient algorithms based on least-squares methodologies. 
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5. Accuracy Performance Between Closed-Formulas and Iterative Procedures in the WSN 
Localization Problem 

In this section we analyze the accuracy performance of both methodologies, optimized and 
closed-formula schemes. Even though the strength of a closed-formula for solving the WSN localization 
problem is its low complexity compared with an iterative algorithm, closed-formulas can present large 
errors in the presence of inaccurate ranging measurements. However, in many cases it is desirable to 
sacrifice accuracy to save energy (i.e., increase battery lifetime). On the other hand, the weakness for 
closed-form methods (i.e., noise sensitivity) represents the strong point for iterative methods and vice 
versa. The goal of both methodologies seems to be in opposite directions. However, the main effort 
in WSN localization research is focused on developing an strategy that can join the strength of both 
methodologies to create an efficient algorithm that can save energy providing the best accuracy in the 
estimated positions. 

Next we present an evaluation of accuracy between closed-formulas and iterative methodologies. 
For the former methods we are considering the classical LS Multilateration, the Min-max method 
(The Min-max approach is based on the intersection of rectangles instead of circle intersections. It 
provides a more simple technique than lateration schemes to obtain position estimates at the expenses 
of accuracy) [38,39], and the bilateration algorithm. For iterative methodologies we are considering two 
algorithms to solve the NLLS: the LM and the TRR algorithms. 

For the simulations that follow, we consider 20 different sensor networks where each one is composed 
by N = 96 unknown sensors, randomly distributed over 100 m by 100 m area. Also, we add M = 4 
non-collinear anchors with full-connectivity on every realization as shown in Figure 4. 

Figure 4. A typical WSN composed by 96 unknown sensors with true positions='<)', 
4 non-collinear anchors ='A', and 96 initial estimates='+'. Each unknown sensor, using 
an initialization algorithm, estimates its initial position by using four reference positions 
(Ai,A2,A3, and A4) and their respective estimated distances. 
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For each network, we add noise to the true distances between anchors and nodes using the log-distance 
path loss model described in Equation (4). The estimated distances are simulated using Osh — 6dB and 
r\p = 2.6, typical parameters for the propagation models on outdoors scenarios, and Pq = —52^ is 
selected according to current commercial specifications for wireless motes [40]. Finally, we assume that 
all nodes have the sensitivity to detect any RF signal coming from anchors. 

To compare the accuracy performance between both methodologies it was necessary to use the same 
set of range measurements for each closed-form method and iterative algorithm. Figure 5 summarizes 
the initial estimates obtained by both methodologies using the RMSE metric as shown the next equation: 

RMSE = J-jy i \ V ?-z i \\ (32) 

where p° represents an initial position estimate for a sensor Si and z ; its true position. 

Figure 5. Position estimates provided by different initialization algorithms. 
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As can be seen, the closed-form LS approach provides the least accurate initial estimates 
(mean = 22.7 m and standard deviation = 2.22 m) compared with iterative algorithms as expected 
due to the noisy ranging measurements. In a similar way, the Min-max scheme also provides large 
errors (mean = 19.7 m and standard deviation = 3.36 m). On the other hand, we can appreciate 
that both iterative algorithms, the LM and the TRR, provide practically the best and similar results 
for initial estimates (mean = 12.54 m and standard deviation = 0.69 m) as expected, and finally the 
bilateration algorithm presents very acceptable initial estimates compared with the last two algorithms 
(mean = 12.96 m and standard deviation = 0.84 m). However, we should consider that the computational 
complexity for the LM and the TRR algorithms is significantly larger than the bilateration algorithm. 
This discussion will be expanded in the next section. 




■ LS (Mean-RMSE= 22.7m, Std=2.22m) 

v Min-max(Mean-RMSE= 19.7m, Std=3.36m) 

• Bilateration (Mean-RMSE=1 2.96m, Std=0.84m) - 

a LM (Mean-RMSE= 12.54m, Std=0.69m) 

> TRR (Mean-RMSE= 12.53m, Std=0.68m) 
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Also, we tested the SX, SI and GSLS algorithms [25] using the same set of networks. The estimated 
positions presented large errors under this scheme as indicated by [26]. Then, these results were 
disregarded in our analysis. 

6. Computational Complexity Analysis between the Bilateration and the LM Algorithm 

The efficiency of an algorithm can be described in terms of the time or space complexity [41]. 
Time complexity refers to the relation between the number of required computations for a given input. 
The space used for a given input provides the space complexity of an algorithm. The computational 
complexity of an algorithm could be described as the number of operations that it takes to find a 
solution [42]. 

In this section we provide an operation count on the number of additions (ADDs), multipliers (MULs), 
divides (DIVs), and square roots (SQRTs) exactly in the way that DSP algorithms are described [43]. 
This will allow an "apple-to-apple" comparison. Moreover, an accurate description lends itself to a cycle 
accurate description for any microprocessor and more significantly, the use of energy models based on 
computing cycles to estimate the energy consumption for a given algorithm. An energy analysis will be 
a discussion of a future work. In next subsections we present the computational complexity analysis for 
the iterative LS and the bilateration algorithm. 

6.1. Computational Analysis of the LM Algorithm 

The LM algorithm could be considered as too expensive for motes given its iterative nature and the 
need to estimate first and second order information {i.e., gradients, Jacobians, and Hessians). The number 
of iterations K is highly dependent on the initial point xo and could be considered a random variable. On 
the other hand, if a good initial estimate, xq, is provided, then the number of iterations is expected to be 
low given the convergence properties of LM. 

We are interested in providing an algorithmic analysis that provides a detailed description in terms 
of additions and subtractions (jointly referred as ADDs), multiplications (MULs), divisions (DIVs), and 
square roots (SQRTs). For simplicity in the next paragraphs we let J^. = J(x k ) and R k = R(xfc). 

The square root is a relevant operation as the error function and the Jacobian estimate requires 
£2 norms to compute distances between sensor and anchors. We also note that the complexity of the 
operations is not the same in terms of the processing resources (hardware and software) they take. 
Abusing notation we have 

ADD < MUL < DIV < SQRT (33) 

This analysis also focuses on the most efficient implementation in terms of the proper operation 
sequencing in order to favor reuse of terms (i.e., avoid computing the same quantity twice). 

We perform the analysis for a single iteration of the LM algorithm, and the total cost for each operation 
is multiplied by K. We also note that K can be modeled as a random variable; the usefulness of this 
approach is discussed later. We assume there are M anchors which have broadcast their position to all 
the nodes. Each node will run the LM algorithm to find its initial position as described in Subsection 3.2. 
We identify three core operations: £2 or Euclidean norm, the error vector and an estimate of J k . 
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The £2 norm will be used to compute the magnitude of the difference between two vectors a, b G I 2 
given by ||a — b|| = a/ (a x — b x ) 2 + (a y — b y ) 2 . This requires three ADDS, two MULs and one SQRT 
The norm is used to compute given in Equation (18) and to estimate the Jacobian as follows: 



-(yi-ofc) 



V (xi-xk) 2 +(yi-yk) 2 V (-xi — -*-fc) 2 +(>'i — y*:) 2 



. V (^M-^) 2 +(>'M->'i : ) 2 V {xni-xk) 2 +(yM-yic) 2 



(34) 



J Mx2 



For Ryt we see that we require M ADDs and M ^2-norms. Accounting for the norms, the error function 
requires AM ADDs, 2M MULs, and M SQRTs. These numbers are recorder in Table 1 . A similar analysis 
follows for Jk- A direct look at Equation (34) indicates that we have the same norm across rows, so we 
can compute them first and then we would need an additional 2M ADDs and 2M DIVs. However, a 
better approach would be to compute the terms l/||xj — x^|| first so that we would require M DIVs, 2M 
MULs, and 2M ADDs. We exchange M DIVs by 2M MULs under the typical case that MULs have a 
much lower complexity than DIVs, particularly for the case of floating point operations. The complexity 
for the Jacobian estimate is also shown on Table 1 . 



Table 1. LM Cost Functions. 



Title 


ADD 


MUL 


DIV 


SQRT 




AM 


2M 


0 


M 


h 


5M 


AM 


M 


M 




3M-3 


3M 


0 


0 


M' f 


2M-2 


2M 


0 


0 


M f 


M-1 


M+l 


0 


0 




3 


3 


0 


1 




3 


6 


1 


0 




2 


4 


0 


0 


Sufficient Decrease 


T(M + A) 


T(M + 2) 


0 


0 


Update 


2 


2 


0 


2 


Stopping Condition 


3 


2 


0 


1 


Total 


(M + 4)T+15M + 7 


(M + 2)r + 12M+18 


M+l 


2M + 4 



Once these two quantities have been evaluated, their use trickles down through the algorithm. The 
costs for the different steps or operations is presented in the remaining part of Table 1 . We just make two 
more remarks on the algorithm complexity. First, note that the approximation to the Hessian matrix J^J^ 



V 2 /(x fc )=j[j, 



M 

I 

.7=1 



( x j- x kf 



({xj-x k ) 2 +(yk-yj) 2 ) 

y I (xj-x k )(yj-y k ) 
7=1 \ ((xj-xk) 2 +(yj-y k )) 




{xj-xk){yj-y k ) 
((xj-xk) 2 +(yj-yk) 2 ) 

(yj-yk) 2 
((xj-xk) 2 +(yk-yj) 2 ) 



(35) 



2x2 
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is of size 2x2 which makes its inversion trivial when computing the LM step Aim, as shown in 
Algorithm 1. 

A LM =(j k J T k +p k iy l J T k R k (36) 
where the gradient of the function J^R/t is given by 

g I ( x j- x k)-(Rkj-\/(xj-x k ) 2 +(yj-y k ) 2 ) \ 
g / (yj-y*)- (^-\/(^-^) 2 +(>7-») 2 ) \ 

j= l \ V ( Jc y- JC *) 2 +(yj-w) 2 / . 2x1 

Second, satisfying the sufficient decrease condition is also an iterative procedure where different 
values of oc^ are tested. We identify T as the number of iterations needed to satisfy this condition. As we 
discuss later, we will model T as a random variable. 

The last row of the table provides the total which we identify as Tadd, Tmul, Tdiv and Tsqrt 
respectively. These numbers are the operations for a single iteration of the LM algorithm. Then, 
for K iterations we have the total number of operations to be Kadd = K ■ Tadd, K M ul = K ■ T MUL , 
Kdiv = K ■ Tdiv, and Ksqrt = K ■ Tsqrt- 

Since the values of T and K are random variables, then a more convenient approach to quantify the 
number of operations would be to look at the average number of operations, i.e., the expected value. It 
is intuitive to assume that T and K are independent, and that for a given network their distributions will 
be identical. Hence, we define 



Kadd = £{K A dd} = e{K}e{T ADD } = e{K}(e{T}(M + 4) + 15M + 7) (38) 

Kmul = £{K M ul} = £{K}£{T MUL } = e{K}(e{T}(M + 2) + 10M+ 15) (39) 

K Drv = e{K DIV } =e{K}(M+l) (40) 

Ksqrt = £{K S qr T } = £{K}(2M + 2) (41) 



where e{x} represents the expected value of the random variable x. Finally, we can quantify the total 
complexity of the LM algorithm by converting operations to a common denominator and compute a 
single representative number that can be used for comparison with other algorithms. The typical way to 
quantify operations is to use the number of processor cycles (on the average) required to complete each 
type of operation. Let us define Nadd, Nmul, Ndiv, and Nsqrt as the number of cycles required for 
floating addition (or subtraction), a multiplication, a division, and square root, respectively. We should 
note that these numbers depend on the micro-processor and hardware used by the mote and the compiler 
tools used to develop the software. Hence, in practice the best way to obtain these values is through code 
profiling using a cycle-accurate simulator. Moreover, as discussed in [44], the number of task cycles can 
be used as part of models that measure energy consumption. Hence, as final measure of complexity for 
the LM algorithm we compute the total number of cycles as 

Nlm = NaddKadd + NmulKmul + Ndiv Kdiv + Nsqrt Ksqrt (42) 



V/(*0 = J[R, = 
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6.2. Computational Analysis of the Bilate ration Algorithm 

The bilateration algorithm is very simple and non-iterative. For M anchors, a sensor node picks (^f) 
pairs of sensors and computes the intersections of the imaginary circles around each anchor with a radius 
given between the anchor and the sensor node. These intersections are computed using geometry with 
a procedure described by Equations (24-28). Then, a cluster with half of the computed intersections is 
found, providing an indication of the area where the node position is located. The number of operations 
required to compute two intersections is presented in Table 2. 



Table 2. Bilateration Cost Operations. 
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Total (2 circle intersections) 


11 


12 


3 


2 


Total Q node combinations 




122 


32 


2Q 



Since this process is repeated Q= ( 2 ) times, then the final row reflects the total operations multiplied 
by this factor. As the intersections are computed, the search for the cluster is performed by Equations (29) 
and (30). Since there are 2Q intersections, we need to select the Q that cluster together (i.e., eliminate 
mirrors). The clustering is based on looking at the distance between all possible pairs of intersections 
and selecting those that exhibit the closets distances among themselves. This requires the calculation 
of S = 2 £( 2 £ -1 ) squared norms, and the use of a clustering or sorting algorithm to find the smallest Q 
elements from the list of S norm values. Taking advantage of the structure of the location points (i.e., the 
two intersections from the same anchor pair are not compared), we can expect an average complexity 
of 0(S) sorting steps using an algorithm like Quickselect algorithm [45]. Hence the final computational 
cost for the bilateration algorithm is presented in Table 3. 



Table 3. Final Computational Cost for the Bilateration Scheme. 



Action 


ADD 


MUL 


DIV 


SQRT 


SORT 


Circle Intersections 


ne 


122 


32 


22 
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Squared Norms 


3S 


2S 


0 


0 


0 


Number of Comparisons 


0 


0 


0 


0 


O(S) 



As with the LM algorithm, we close this subsection by providing an expression in terms of processor 
cycles. Using the same characterization for all main operations of the algorithm, we can provide a total 
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cycle count that can be directly compared with other algorithms. Obviously, a lower cycle implies lower 
complexity when the hardware and software development tools are identical. The expression for total 
cycles is 

N BL = N ADD ■ (1 1 Q + 35) + Nmul • ( 1 2Q + 25) + N DIV ■ (3 Q) + N S qr T -(2Q) + N SO rt (43) 

It is easy to see that the bilateration scheme uses a significantly less number of cycle for all operations. 
Experimental data in [45] indicates that the cycle count for the complete sorting step with the QuickSelect 
algorithm with a pipelined architecture can be achieved within 2500 and 3000 clock cycles. 

To complete the computational complexity analysis, we need the number of CPU cycles required 
for the four basic operations as floating point operations. These values are highly dependent on the 
architecture of the mote processor. An extensive study in [46] provides good representative values for 
processors with some level of hardware support. The values are summarized on Table 4, and it shows 
the relation between basic operations and CPU cycles. 



Table 4. Operation cycle counts. 



ADD 


MUL 


DIV 


SQRT 


11 


25 


112 


119 



Next, we use Tables 1 and 2 to obtain the number of CPU cycles required by each initialization stage, 
BL and LM respectively. For the LM initialization stage we are using M = 4 anchors and the random 
variables T and k. We recall that k is the number of iterations spent by the LM algorithm to find a 
solution. These values are obtained through simulations where e{T} = 2 and e{k} = 13. In this way the 
total cycles required by the LM algorithm according to Equation (42) is given by 

Kim =e{£}[((M+4)e{r} + 15M + 7) (11)+ 
((M + 2)e{r} + 12M+18)(25)+ 

(M+l)(112)+ (U) 
(2M+ 4) (1 19) ] = 63063 cycles 

Similarly, the total number of cycles used by the BL stage is given by Equation (43) as 

K BL =(U)(UQ + 3S) + 
(25)(12g + 25) + 

(112)(3g)+ (45) 

(119)(2g)+ 

(Nsort) = 141 98 cycles 

where 2 = 6,5 = 66, and Nsort = 2750. The value for Nsort represents the total number of cycles 
required to perform the sorting step of the BL algorithm. This step can be performed using efficiently the 
Quickselect algorithm [47]. As expected, the LM algorithm consumes more energy in the initialization 



Sensors 2012, 12 



859 



process than the BL scheme. However, the former represents a better choice when accuracy is required. 
Therefore, the BL can be an alternative localization scheme when a good tradeoff between accuracy and 
energy consumption is required on the initial estimates. 

7. Conclusions 

In this research, we analyzed a localization algorithm that can be realistically deployed over real 
WSNs which can provide good accuracy performance with low computational complexity. The 
bilateration algorithm is a distributed scheme that can be used as an initialization stage to find an initial 
set of locations. 

Most initialization algorithms demand very high computing power to provide a set of initial estimates 
for an N-node WSN. The analyzed algorithm is capable to provide competitive initial estimates at 
low processing power. This approach is basically formed by two stages. The first stage consists of 
finding all circle intersections formed by anchor positions and their respective range estimates to a 
sensor node, obtained by ranging techniques like ToA or RSS. The great advantage of this approach 
is to use "closed-formulas" to find all circle intersections (i.e., candidate positions) using two anchors 
at a time. In the second stage, the algorithm uses a sorting algorithm to find the cluster of candidate 
positions that tend to be closer to each other around the true location. The cluster with the nearby 
candidate positions is averaged to finally obtain the initial location. This scheme can be used by any 
WSN localization algorithm that needs initial approximations. Also, it is implementable in constrained 
devices with low processing and memory capabilities (i.e., motes). Results show that this initialization 
algorithm is well behaved (e.g., computational and accuracy performance) in comparison with other well 
known algorithms like LS methodologies. 

Finally, we are interested in exploring iteratively, at the refinement process, the Levenberg-Marquardt 
approach for node localization. We believe that this methodology can play a crucial role in producing 
excellent position estimates with high accuracy and low energy consumption due to the rate of 
convergence associated with this optimization technique. 
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